function [min_male, max_male, min_female, max_female, problem] = get_ceb(couple,singlemale,singlefemale,leisurem,leisuref,hworkm,hworkf,wagem,wagef,q,Q,nonlabor,mutually_acceptable,divorce_costsM,divorce_costsF,divorce_costsMF,text_solver)

% ================
% define constants
% ================

T = 112;                                              % total potential working hours
H = length(leisurem);                                 % total number of households in the data
alpha = 0.4;                                          % bound for non-labor income

% =========================
% define decision variables
% =========================
yalmip('clear')
qM = sdpvar(H,1);                                     % private consumption of Hicksian good male
qF = sdpvar(H,1);                                     % private consumption of Hicksian good female

pM = sdpvar(H,H,'full');                              % personalized public good price male
pF = sdpvar(H,H,'full');                              % personalized public good price female

pmhworkM = sdpvar(H,H,'full');                        % personalized male's hwork price male
pmhworkF = sdpvar(H,H,'full');                        % personalized male's hwork price female

pfhworkM = sdpvar(H,H,'full');                        % personalized female's hwork price male
pfhworkF = sdpvar(H,H,'full');                        % personalized female's hwork price female

wageM = sdpvar(H,1);                                  % shadow wage male
wageF = sdpvar(H,1);                                  % shadow wage female

nl_incomeM = sdpvar(H,1);                             % non-labor income male
nl_incomeF = sdpvar(H,1);                             % non-labor income female

sharem = sdpvar(H,1);                                 % consumption based RICEB male
sharef = sdpvar(H,1);                                 % consumption based RICEB female


%%
% ==================
% define constraints
% ==================

Constraints = [];

% restriction on individual private consumption
for i=1:H    
    Constraints = [Constraints, q(i) == qM(i) + qF(i), qM(i) >=0, qF(i) >= 0];
    
    if singlefemale(i) == 1
        Constraints = [Constraints,qM(i) == 0];
    elseif singlemale(i) == 1
        Constraints = [Constraints,qF(i) == 0];
    end
    
end


% shadow price should be non-negative and equal to observed wage if the individual is employed
for i=1:H
    
    Constraints = [Constraints,wageF(i) >= 0, wageM(i) >= 0];
    
    if ~isnan(wagem(i))
        Constraints = [Constraints, wageM(i) == wagem(i)];
    end
    if ~isnan(wagef(i))
        Constraints = [Constraints, wageF(i) == wagef(i)];
    end
    
end


% restriction on public prices
for i=1:H
    for j=1:H
        
        Constraints = [Constraints, 1 == pM(i,j) + pF(i,j), pM(i,j) >=0, pF(i,j) >= 0];
        Constraints = [Constraints, wageM(i) == pmhworkM(i,j) + pmhworkF(i,j), pmhworkM(i,j) >=0, pmhworkF(i,j) >= 0];
        Constraints = [Constraints, wageF(j) == pfhworkM(i,j) + pfhworkF(i,j), pfhworkM(i,j) >=0, pfhworkF(i,j) >= 0];
        
        if (couple(i) == 1 || singlemale(i) == 1) && singlemale(j) == 1
            Constraints = [Constraints, 1 == pM(i,j), wageM(i) == pmhworkM(i,j), wageF(j) == pfhworkM(i,j)];
        end
        
        if (couple(j) == 1 || singlefemale(j) == 1) && singlefemale(i) == 1
            Constraints = [Constraints, 1 == pF(i,j), wageM(i) == pmhworkF(i,j), wageF(j) == pfhworkF(i,j) ];
        end
    end
    
end


% feasibility restrictions on the individual nonlabor income (for each matched pair)
% individual nonlabor incomes after divorce must lie between 40% and 60% of total nonlabor income under marriage
for i=1:H
    
    Constraints = [Constraints, nonlabor(i) == nl_incomeM(i) + nl_incomeF(i)];
    
    if couple(i) == 1
        if nonlabor(i)>=0
            Constraints = [Constraints, alpha*nonlabor(i) <= nl_incomeM(i), nl_incomeM(i) <= (1-alpha)*nonlabor(i), alpha*nonlabor(i) <= nl_incomeF(i), nl_incomeF(i) <= (1-alpha)*nonlabor(i)];
        else
            Constraints = [Constraints,(1-alpha)*nonlabor(i) <= nl_incomeM(i), nl_incomeM(i) <= alpha*nonlabor(i), (1-alpha)*nonlabor(i) <= nl_incomeF(i), nl_incomeF(i) <= alpha*nonlabor(i)];
        end
    elseif singlefemale(i) == 1
        Constraints = [Constraints, nonlabor(i) == nl_incomeF(i)];
    elseif singlemale(i) == 1
        Constraints = [Constraints, nonlabor(i) == nl_incomeM(i)];
    end
    
end


% individual rationality 
for i = 1:H
        if couple(i) == 1        
            Constraints = [Constraints, wageM(i)*T + nl_incomeM(i) - divorce_costsM(i) <= qM(i) + wageM(i)*leisurem(i) + Q(i) + wageM(i)*hworkm(i) + wageF(i)*hworkf(i)];
            Constraints = [Constraints, wageF(i)*T + nl_incomeF(i) - divorce_costsF(i) <= qF(i) + wageF(i)*leisuref(i) + Q(i) + wageM(i)*hworkm(i) + wageF(i)*hworkf(i)];        
        end        
end
    

% no blocking pairs
for i = 1:H
    for j = 1:H        
        if mutually_acceptable(i,j) == 1
            Constraints = [Constraints, (wageM(i) + wageF(j))*T + nl_incomeM(i) + nl_incomeF(j) - divorce_costsMF(i,j) <= ...
                qM(i) + qF(j) + wageM(i)*leisurem(i) + wageF(j)*leisuref(j) + pM(i,j)*Q(i) + pF(i,j)*Q(j) + ...
                pmhworkM(i,j)*hworkm(i) + pmhworkF(i,j)*hworkm(j) + pfhworkM(i,j)*hworkf(i) + pfhworkF(i,j)*hworkf(j)];
        end
        
    end
end


% define RICEBs
for i = 1:H
    if couple(i) == 1
        Constraints = [Constraints, sharem(i) == qM(i) + Q(i),...
            sharef(i) == qF(i) + Q(i)];
    elseif singlemale(i) == 1
        Constraints = [Constraints, sharem(i) == qM(i) + Q(i), sharef(i) == 0];
    elseif singlefemale(i) == 1
        Constraints = [Constraints, sharem(i) == 0, sharef(i) == qF(i) + Q(i)];
    end
end



%%
% =================
% Solve the problem
% =================

min_male = zeros(H,1);
max_male = zeros(H,1);
min_female = zeros(H,1);
max_female = zeros(H,1);
problem = 0;

param = eye(H,H);
u = sdpvar(H,1);    % vector of parameters necessasry to specify the function
options = sdpsettings('verbose',0,'solver',text_solver);

Objective1 = sum(u.*sharem);
Opt1 = optimizer(Constraints,Objective1,options,u,sharem);

Objective2 = sum(u.*sharef);
Opt2 = optimizer(Constraints,Objective2,options,u,sharef);

for i=1:H
    
    [xvalue, errorcode] = Opt1(param(:,i));
    if errorcode~=0
        min_male(i) = NaN;
        problem = 1;
    else
        min_male(i) = xvalue(i);
    end
    
    [xvalue, errorcode] = Opt1(-param(:,i));
    if errorcode~=0
        max_male(i) = NaN;
        problem = 1;
    else
        max_male(i) = xvalue(i);
    end

    
    [xvalue, errorcode] = Opt2(param(:,i));
    if errorcode~=0
        min_female(i) = NaN;
        problem = 1;
    else
        min_female(i) = xvalue(i);
    end

    
    [xvalue, errorcode] = Opt2(-param(:,i));
    if errorcode~=0
        max_sharef(i) = NaN;
        problem = 1;
    else
        max_female(i) = xvalue(i);
    end

end


return


